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The method of phenomenological damping developed by Pitaevskii for superfluidity near the A 
point is simulated numerically for the case of a dilute, alkali, inhomogeneous Bose-condensed gas 
near absolute zero. We study several features of this method in describing the damping of excitations 
in a Bose-Einstein condensate. In addition, we show that the method may be employed to obtain 
numerically accurate ground states for a variety of trap potentials. 
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Recent experiments with dilute, trapped Bose-Einstein condensates (BEC) have indicated the presence of damping 
of collective excitations The most commonly used theoretical tool so far for describing BEC, the Gross-Pitaevskii 
equation (GPE) Q|, has been shown to give accurate predictions for the frequencies of elementary excitations ||-[5), 
but does not contain any mechanism for the damping of these excitations. Theories that go beyond the GPE are now 
subjects of intense study j7u|. One of the theories which incorporates dissipation in a quantum fluid is Pitaevskii's 
phenomenological theory of superfluidity near the A point ||], and we have used it to investigate damping in a dilute, 
trapped BEC numerically. This phenomenological method of including damping should, as argued by Pitaevskii, be 
valid when the system is close to thermal equilibrium. We shall begin by deriving the required modification of GPE. 
We intend to demonstrate in this paper two distinct uses of the equation: its use in describing damping of oscillations 
in trapped atomic BEC, and its use as a computational method to generate accurate cigcnstates of the standard GPE. 
—i ■ The standard GPE describes the motion of the mean field, if), which is taken to be the ensemble average of the 
' boson field operator We write the GPE in the form 
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, 1 1 , where N is the number of particles, a the interatomic s-wave scattering length and m the mass of a single particle. 
The mean field tp is customarily taken to be the condensate wave function, although it is also consistent to define ip 
to be the condensate ground state plus any coherent excitation. 

From inspection, the GPE contains no term which can describe damping. The first two terms are simply the 
energy of a single particle in the trap while the third describes the non-dissipative effect of other particles in the 
system. Attempts to improve on the GPE include taking into account the correlation effects of excitations and the 
introduction of the many-body T-matrix [g . These theories so far have not produced any explicit results for damping 
in the inhomogeneous case, although recent related work by Giorgini has shown the presence of Landau and Beliaev 
damping pfj[ . 

It is, however, possible to include damping phenomenologically, in a manner similar to that introduced by Pitaevskii 
H . The argument which we present below to obtain our expression for phenomenological damping follows this work 
closely. The procedure described is, in fact, a widely used one in the more general problems of dynamical critical 
phenomena Q. 

First of all, we note that any damping process eventually leads to an equilibrium state. For the condensate at T = 
such an equilibrium is described by the equation 

h 2 

~ + VtrapW^ + C„HV ~^ = 0. (3) 

This expression may be derived from minimizing an appropriate energy functional or equivalently by finding the 
eigensolution to the time dependent GPE. Because the process we want to describe is a relaxation process, in our 
equation of motion for the condensate ip, 

ih^ = £V, (4) 
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the operator C cannot be Hermitian. The anti-Hermitian part of C is associated with the processes by which equilib- 
rium is approached. This anti-Hermitian part must vanish at equilibrium, that is, when Eq. (0) is satisfied. For the 
cases where the deviation from equilibrium is small, one may therefore write the anti-Hermitian part in the form 

^A Q^VSV + ^trap(r)^ + C n |V|V ~ > ( 5 ) 

where A is a dimensionless factor which is inversely proportional to the relaxation time. The Hermitian part is given 
by the fact that for A = 0, the theory must reduce to the conventional GPE at T = 0; the final equation including 
relaxation to equilibrium is, then, given by 

= ( 1 + iA ) (~^- V ? + WO + C„|Vf - im) V, (6) 



dt |_ 2m 

where A < for damping. 

Equation (^) is, in fact, a rather general equation of motion which describes evolution towards equilibrium. Damping 
of elementary excitations may be described in the frame work of Eq. (^|) by considering our mean field ip to be broken 
into two parts: the condensate ground state tp g and a coherent excitation, 5, 

^ = e- i "*(^ + 5). (7) 

Evolution in Eq. @ results in the damping of 6 towards zero. The exact value for the damping parameter A is 
expected to depend on a variety of factors. One would, for instance, expect it to be temperature dependent, since the 
thermal component present in the trap is clearly going to be one of the main sources of damping. In the original case 
of absorption of sound in supcrfluid Helium considered by Pitaevskii, A could be expressed in terms of the coefficients 
of second viscosity used by Khalatnikov jl2| . A closely related result (which actually derives the Landau-Khalatnikov 
two-fluid model for the homogeneous Bose-condensed gas) has been calculated by Kirkpatrick and Dorfman using the 
Chapman- Enskog procedure |i~3| ] . They also explicitly derived various transport coefficients such as the viscosity for 
their homogeneous case. At present, a complete analysis for the inhomogeneous case is not yet available. 

Some rough initial estimate for the damping coefficient A may, however, be made by noting that, physically, A 
represents the rate at which the excited components turn into the condensate. This may be approximated by the 



transition probability W + (N) that has been estimated in a recent work by Gardiner et al. 14 using the quantum 
kinetic theory. W + (N) gives the rate at which the thermal particles above the condensate band enter the condensate 
due to collisions. Although the description of the process of evaporative cooling was their primary goal, their main 
assumption in formulating W + (N) was that the condensate does not readily act back on the thermal component 
to change its temperature. This assumption is expected to be valid when one considers quasiequilibrium situations. 
Using the reported values from the MIT experiment || as an example, and the expression 

where K\{z) is a modified Bessel function, we obtain A ~ W + {N)/lo « —0.03, for T « T c /10 and a ss 3.45nm. Here, 
w was taken to be the trap frequency relevant to the oscillatory motion, that is u = 2ir x 19Hz. The corresponding 
damping time is then approximately 200 ms, which is of the same order of magnitude as the experimental damping 
time of 250 ms [§. 

Numerical simulation of Eq. is done using a 4th order Runge-Kutta algorithm. Although the simulations in 
this paper are performed using a one dimensional (ID) Runge-Kutta integrator, the general arguments that follow do 
not depend on dimensionality, and it is straight forward to apply the ideas of the present paper to a 3D simulation. 
We simulate damping of elementary excitations by first populating a mode of elementary excitation on a condensate 
and then time-evolving this condensate plus excitation using Eq. (^|). In practice, such a state is generated by, for 
instance, applying time dependent perturbations to the confining potential initially and then observing the subsequent 
behaviour of the condensate after turning off the perturbation. Exactly how the excitations damp using this formalism 
can be analysed quantitatively using the method of quasiparticle projection |l5| ]. This projection method is based on 
the Bogoliubov transformation of fluctuations about the condensate mean field: one writes a general wave function tp 
as 
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where the index i denotes quasiparticle energy levels, and the functions Uj(r) and Uj(r) obey a set of coupled equations 
usually known as the Bogoliubov-de Gennes equations which diagonalise a quadratic Hamiltonian |l6f| . The sum over 
i on the right hand side of Eq. (||) represents the fluctuation about the condensate ground state ip g which is itself 
weighted by a factor (1 + b g ). The population of quasiparticles in level i is then given by the quantity b*bi = |6;| 2 
where the coefficients hi are taken to be the mean values of the quasiparticle annihilation operators. By using the 
well known orthogonality relations of Ui(r) and i>i(r), one may extract the expansion coefficients b „ a nd bi uniquely 
at any point in time, and hence deduce the evolution of the quasiparticle population in each mode |lq ]. 

In order to simulate the experimental results as closely as possible, we initially populate the breathing mode, the 
second lowest excitation, by an arbitrary amount and then plot the evolution of this population, |i>2| 2 , over time. 
Figure 1 shows the population of mode 2 over time as well as the population in the condensate mode \b g \ 2 for various 
values of A. An initial population of 0.25 means | £> 2 1 2 = 0.25. The corresponding ground state population is shown in 
the lower part of the plot. We see that for smaller values of A, i.e. when the excited component does not damp out 
quickly over time, there is a greater interplay between the ground and the excited state populations due to nonlinear 
mixing. 

The width of the condensate is indeed what is actually measured in the experiment, and therefore this is plotted as 
a function of time in Fig. 2. The general shape reported in Ref. of the form Aexp(— Ai)sin(27wi + </>) + B is clearly 
reflected in this figure. The solid curve gives what is expected for the MIT experiment [gj (i.e. A w —0.03). As the 
simulation was done in ID, however, we do not intend to claim any quantitative analogy. Also, in real experiments, 
several quasiparticle modes, rather than just the breathing mode are simultaneously excited. It is clear, however, 
that this description of damping should be satisfactory in cases where an exact formalism is not needed: for instance, 
when A may be treated as a parameter to be obtained from an experiment. 

One of the ways to characterise different descriptions of damping is to see its dependence on energy. For instance, 
in the homogeneous gas, it is well known that the damping rate increases as k 2 when the relaxation process is due 
to thermal conduction and viscosity pd[ |. We note that the Pitaevskii prescription gives a damping rate which has 
a linear dependence on the mode energy, with the high energy components being damped out more rapidly. This 
aspect may be seen simply by inspection of the equation of motion Eq. (^), and it can also be confirmed by explicitly 
calculating the damping rate of the quasiparticle populations, |6i| 2 , for various values of i by linearisation of Eq. (|^). 

We now look at the second possible use of Eq. (g) . We note that finding a reliable solution to the time- independent 
GPE is an essential first step in many numerical simulations involving a trapped BEC. The Pitaevskii damping 
description may, in fact, be used to this end as an effective numerical tool to find the eigenstate. This is possible 
because the equilibrium state Eq. (|^) is itself a solution to Eq. (|^). To find the eigenstate, one starts by initially 
choosing an arbitrary function (preferably close to the desired solution for faster convergence, e.g. for a condensate 
in a trap, a Gaussian or a Thomas-Fermi solution) and run the routine with pre-determined values of /i and A. It 
should be noted here that, since the equation is non-hermitian, the evolution does not conserve the norm of ip. It is 
therefore necessary to proceed by first running the routine with an estimate of C n for the chosen fj, until convergence, 
and then later adjusting C n from the renormalisation of tp, 

J \4>(v)\ 2 dv=l. (10) 

The correct C n for the chosen is then given by the initial estimate divided by the final normalisation constant. 

We found that the time it takes to find a solution in this way is very much shorter than other methods such as the 
adiabatic turn on of the non-linearity C n . As a rough indication, evolution time of t = IQ/uj was required for growing 
an eigenstate with C n = 100 and A = — 2 when using the damping method while it took us t = 500/a; using the 
adiabatic method; in real time on a Sun SPARC station 10, these running times corresponded to roughly 10 minutes 
and 3 hours respectively. The damping method is able to find a solution for as large a C n as desired, provided there 
are a sufficient number of grid points to accommodate the larger condensate. An exact eigensolution to GPE has a 
flat phase profile; the damping simulation with A = —2, C n = 50 running for t = 15/uj typically gave a maximum 
change in phase across the condensate of the order of 10~ 7 radians, indicating the high accuracy of this method. 

Other advantages of this approach include its flexibility in handling traps of any arbitrary shape, and the fact 
that it is not necessary to expand the wave function in a set of complete, orthogonal basis states which can be a 
computationally expensive operation. 

The value of A cannot, however, be increased indefinitely for faster convergence, as this results in numerical 
instability. It was found that we get numerical instability when Ai s t cp > 0.03 where t s t ep is the size of the incremental 
time step in numerical integration. The onset of numerical instability is independent of the total length of simulation 
and the size of the non-linearity constant C n . One may understand this from the fact that the damping per time step 
is roughly given by e _Atatop so that e~ ' 03 w 0.97 implies a reduction in the amplitude by about 3 % per numerical 
time step. Since the typical step size used is of the order of 10~ 3 harmonic oscillator units, this, in fact, represents a 
complete decay on a time scale of order a period. 
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We also point out that the vortex eigenstate in a trap may be obtained numerically by using this damping pre- 
scription; one needs only to start with the first excited harmonic oscillator state which has an odd parity. One of the 
corollaries of the damping method is that although parity is conserved the number of zero crossings is not; if we started 
with, say, the 3rd excited harmonic oscillator state, the resulting solution is still the state with one zero crossing. The 
phase profile of such a condensate displays a clear tt phase jump at the origin, while being flat everywhere else. 

In summary, we found that the phcnomcnological damping description docs indeed give results which are compatible 
with those of the experiments. In a more general context, it provides an efficient numerical tool for finding the 
eigenstate of the time independent GPE. 
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Evolution of quasiparticle and condensate population 
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FIG. 1. Population of the breathing mode, |i>2| 2 , over time for various values of A. Solid, dashed, dot-dashed, and dotted lines indicate 
A of —0.03, —0.1, —0.25, and —0.5 respectively. The lower half of the plot shows the corresponding variation in |l> 9 | 2 . 
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FIG. 2. Variation of the width of the trapped condensate over time. The width is measured by the standard deviation of x, the 
position, weighted by the corresponding condensate density. Solid, dashed, dot-dashed and dotted lines indicate A of —0.03, —0.1, —0.25, 
and —0.5 respectively, as in the previous figure. 
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